Highly accurate and efficient self-force computation using time-domain methods: 
Error estimates, validation, and optimization 
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If a small "particle" of mass ^lM (with ^ <C 1) orbits a Schwarzschild or Kerr black hole of mass 
M, the particle is subject to an 0(/i) radiation-reaction "self-force". Here I argue that it's valuable 
to compute this self-force highly accurately (relative error of < 10"'') and efficiently, and I describe 
techniques for doing this and for obtaining and validating error estimates for the computation. I 
use an adaptive-mesh-refinement (AMR) time-domain numerical integration of the perturbation 
equations in the Barack-Ori mode-sum regularization formalism; this is efficient, yet allows easy 
generalization to arbitrary particle orbits. I focus on the model problem of a scalar particle in a 
circular geodesic orbit in Schwarzschild spacetime. 

The mode-sum formalism gives the self-force as an infinite sum of regularized spherical-harmonic 
modes X^^o-^^'^cg, with -Fi,rog (and an "internal" error estimate) computed numerically for I < 
30 and estimated for larger I by fitting an asymptotic "tail" series. Here I validate the internal 
error estimates for the individual -Fi,reg using a large set of numerical self- force computations of 
widely-varying accuracies. I present numerical evidence that the actual numerical errors in -Fi,rcg 
for different I are at most weakly correlated, so the usual statistical error estimates are valid for 
computing the self-force. I show that the tail fit is numerically ill-conditioned, but this can be 
mostly alleviated by renormalizing the basis functions to have similar magnitudes. 

Using AMR, fixed mesh refinement, and extended- precision fioating-point arithmetic, I obtain the 
(contravariant) radial component of the self-force for a particle in a circular geodesic orbit of areal 
radius r — lOM to within 1 ppm relative error, as estimated both by internal error estimates and 
by comparison with previously-published frequency-domain calculations. 

PACS numbers: 04.25.Nx, 04.25.dg 02.70.-c, 04.25. Dm, 

Keywords: self-force, radiation reaction, extreme mass-ratio inspiral, Barack-Ori mode-sum regularization, 
black holes, least-squares fitting, ill-conditioning 



> 

00 
00 

o 
o 



X 



This paper is dedicated to the memory of 
Thomas Radke, my late friend, colleague, and 
partner in many computational adventures. 



I. INTRODUCTION 

Consider a small "particle" of mass /iA/ (with /i <C 
1) moving freely in an asymptotically- flat background 
spacetime, say for definiteness Schwarzschild or Kerr 
spacetime of mass M. This system emits gravita- 
tional waves (GWs), and there is a corresponding 
radiation-reaction influence on the particle's motion. 
Self-consistently calculating this motion and the emitted 
gravitational radiation is a long-standing research ques- 
tion, and is interesting both as an abstract problem in 
general relativity and as an essential prerequisite for the 
success of the proposed Laser Interferometer Space Array 
(LISA) space-based gravitational radiation detector. A 
typical LISA extreme mass ratio inspiral (EMRI) source 
is expected to comprise a stellar-mass black hole or neu- 
tron star (the "particle") orbiting a supermassive black 
hole with M ~ lO'^ MqQ so that n ~ 10"^ to IQ-^; the 



particle's orbit will typically be both inclined (with re- 
spect to the supermassive black hole's equatorial plane 
defined by its spin) and modcrately-to-highly eccentric. 
LISA is expected to observe many such systems, some of 
them at quite high signal/noise ratios once the raw data 
stream is matched- fi ltered against appropriate waveform 
templates ( 124l4l27 |: see section lTl A ll for further discus- 
sion). 

The particle's orbit may be highly r elativ istic, so post- 
Newtonian methods (see, for example, |l28l section 6.10]; 
129l4l32l | and references therein) may not be accurate for 
this problem. Since the timescale for radiation reaction 
to shrink the orbit is very long (^ ijl~^M) while the re- 
quired resolution near the particle is very high (~ uM), 
full numerical- relativity methods (see, for example, |133l - 
Il37l | and references therein) are prohibitively expensive 
for this problemQ 

Instead, it's appropriate to use black hole perturba- 
tion theory, treating the particle as an 0(/i) perturba- 
tion on the background Schwarzschild or Kerr space- 



* IjthornQastro. i ndiana.edu 

^ Here Mq denotes the solar mass. 



^ A number of researchers have attempted to develop special 
numerical-relativity methods to make such simulations practical, 
at least for systems with "intermediate" mass ratios fi ~ 10~^. 
Although promisin g initial results have been obtained (see, for 
example, Il38l4l4a ). it has not (yet) been possible to perform 
numerical evolutions lasting for radiation-reaction time scales. 
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time. A self-consistent calculation of the emitted gravita- 
tional radiation requires knowledge of the metric pertur- 
bation indu ced by the part icle u p to and inclu ding 0(u^) 



terms ( jl43l section 5.5.6]; jl44 section 11.1] 
The theoretical formalism for such calculations is not yet 
fully developed!! here I present calculations only for the 
0(/i) self force. 

Building on the early work of DcW itt and Brehme [l5l| 
(with a correction by Hobbs [l5^)l^ the 0(/i) "MiSa- 
TaQuWa" equations of motion for a gravitational point 
particle in a (strong-field) curved spacetime were first 
derived by Min o, Sasaki, and Tanaka |l54{ and Quinn 
and Wald jl55l ] (see also Detweiler's analysis jl56| ]). and 
have recently been rederived in a more rigorous mai mer 
by GraUa and Wald ^E^E See ^M, [HI ITHfll - flBH for 
general reviews of the self-force problem. 

The particle's motion may be modelled as ei- 
ther (i) non-geodesic motion in the background 
Schwarzschild/Kerr spacetime under the influence of a 
radiation-reaction "self-force" , or (ii) geodesic motion in 
a pert urbe d spacetime. These two perspectives are equiv- 
alent [lei]; in this work I use (i). The MiSaTaQuWa 
equations then give the self-force in terms of (the gradi- 
ent of) the metric perturbation due to the particle, which 
must be computed using black-hole perturbation theory. 

The computation of the metric perturbation due to a 
point particle is particularly difficult because the "per- 
turbation" is formally infinite at the particle. A practi- 
cal "mode-sum" scheme to regularize the metric pe rtur- 
bation was developed by Barack and Ori |l63l4l67| . and 
in slightl y different forms by Detweiler, Messaritaki, and 
Whiting [m [16^ and Haas and Poisson ^7^. Here I 
follow the B arack -Ori "£-mode" regularization (described 
in detail in jl67| and summarized in section llll|) . This 
is based on a spherical-harmonic decomposition of the 
metric perturbation, allowing the 4-vector self-force F°' 
to be written as an infinite sum of regularized modes 
F° = J2iLo ^e.rcg- Each regularized mode F^^^-cg is calcu- 
lated by solving a set of linear partial differential equa- 
tions (PDEs), computing certain derivatives of the PDE 
solutions along the particle worldline, and finally sub- 
tracting certain analytically-known regularization coeffi- 
cients. 

Depending on how the PDEs arc solved, there are two 
broad classes of self-force computations within the mode- 
sum regularization framework: frequency-domain and 
time-domain. Frequency-domain computations involve a 
Fourier transform of each mode's PDEs in time, reducing 



^ See, for example, Il47l4l5dll for recent work towards 0(^^) cal- 
culations. 

* Another significant early work is that of Gal'tsov |153|| . but this 
approach has serious causality difficulties: in a curved spacetime 
it gives the self-force at a specified time in terms of the future 
evolution of the particle. 

^ Gralla, Harte, and Wald |158| have also recently obtained a rigor- 
ous derivation of the electromagnetic self-force in a curved space- 
time. 



the numerical computation to the solution of a set of or- 
dinary differential equations (ODEs) for each mode (see, 
for example, jl69l ]). Frequency-domain computations are 
typically very efficient and accurate for circular or near- 
circular particle orbitslf] but degrade rapidly in efficiency 
with increasing eccentricity of the particle' s orb i t, be com- 
ing impractical for highly eccentric orbits jl72l Il73l ] 13 In 
contrast, time-domain computations involve a direct nu- 
merical integration of each mode's PDEs, and have tra- 
ditionally been somewhat less efficient and accurate than 
frequency-domain computations. However, time-domain 
computations can accommodate arbitrary particle orbits 
with only minor penalties in performance and accuracy 
( |l75l ] ) , and some c omplicati ons in the numerical schemes 
(see, for example, [iTa [iTtI ] ) . 

In this work I use the time-domain approach, using 
an adaptive mesh refi neme nt (AMR) code with 4th or- 
der finite differencing jl78l ] to solve each mode's PDEs 
very accurately and efficiently. To simplify the boundary 
treatment, I use a characteristic (double-null) evolution 
scheme. I restrict consideration to the model problem 
of computing the self-force on a scalar particle moving in 
Schwarzschild spacetime. This is a widely-used test prob- 
lem in the field of self- force cal culations , with past nu - 
merical computations including [HI [HI [ill ITtqI - 
Il86l ] [f|[l For the numerical computations presented here, 
I further restrict consideration to the computation of the 
radial component of the self-force for a scalar particle in 
a circular geodesic orbit about the Schwarzschild black 
hole. However, I also simulate the accuracy to be ex- 
pected when similar methods arc applied to generic non- 
circular particle orbits. 

The basic mode-sum technique for self-force compu- 
tation discussed here is already well-known. The main 
new results in this paper concern (a) the (small) ex- 
tension of these techniques to accommodate the use of 
characteristic AMR for the numerical integrations, (b) 
the error estimates for such a self-force computation, (c) 
the validation of these error estimates using a large set 
of numerical computations of widely-varying accuracies, 
(d) the tail fit's ill-conditioning, (e) the cost/accuracy 
tradeoffs for the computation, and (f) the demonstra- 
tion of consistency at levels of ~ 0.1 parts per million 



® As a notable example of this accuracy, Blanchet et al. [l7ll | 
have recently computed the gravitational self-force for circular 
geodesic orbits in Schwarzschild spacetime to a relative accuracy 
of approximately one p art in 10^^ . 

^ Barack, Ori, and Sago [174( have recently found an elegant so- 
lution for some other limitations which had previously affected 
frequency-domain calculations. 

* The electromagnetic self-force (a more complicated "toy model" 
by virtue of the nontrivial gauge freedom) has been studied 
by Il87l . (Note also the recent work described in footnote [5] ) 
The gravitation al se l f-force has been studie d by numerous au- 
thors, including [ml . [l73l , [iTTI. [Tsl . [187^191 . 

® Warburton and Barack |192|| have recently reported results for 
the self-force on a scalar charge in a circular equatorial geodesic 
orbit in Kerr spacetime. 
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(ppm) relative error between the time-domain self-force 
computations presented here and the highly-accurate 
frequency-domain computations of Detweilcr, Messari- 
taki, and Whiting \l6^ . 

The remainder of this paper is organized as follows: 
Section [I Al outlines the notation used in this paper. Sec- 
tion [n] discusses the scientific importance of highly ac- 
curate and efficient self-force computations. Section IIIII 
outlines the Barack-Ori mode-sum regularization proce- 
dure for self- force computations. Section ITVl outlines the 
numerical methods I use for the self-force calculation and 
its error estimates. Section |V] presents my numerical re- 
sults. Section IVII presents conclusions and directions for 
further research. 



A. Notation 

I ge nera lly follow the sign and notation conventions of 
Wald |193j . with G ^ c = 1 units and a (— , +, +, +) met- 
ric signature. I use the Penrose abstract-index notation, 
with Latin indices ab running over spacetimc coordinates. 
g is the determinant of the 4-metric and Vq the associ- 
ated covariant derivative operator. □ = VqV" is the 
4-dimensional wave operator. || • Umis is the root-mean- 
square norm on ||{xfc}|| 



II. THE IMPORTANCE OF HIGHLY 
ACCURATE AND EFFICIENT SELF-FORCE 
CALCULATIONS 

In this section I outline several different lines of argu- 
ment suggesting that it's scientifically valuable to com- 
pute the EMRI self-force highly accurately and efficiently. 



A. The Importance of High Accuracy 

1. LISA 

A major part of the motivation for self-force cal- 
culations comes from their planned application to 
EMRI data a naly sis for LISA. In the words of Amaro- 
Seoane et al. |l26l . section 3.1], 

A typical EMRI signal will have an instanta- 
neous amplitude an order of magnitude be- 
low the LISA'S instrumental noise and (at 
low frequencies) as many as several orders of 
magnitude below the gravitational wave fore- 
ground from Galactic compact binaries. This 
makes detection a rather difficult problem. 
However, the signals are very long lived, and 
will be observed over more than 10^ cycles, 
which in principle allows the signal-to-noise 
ratio (SNR) to be built up over time using 
matched filtering. 



Matched filtering of the entire years-long LISA data 
stream would be impractically expensive for det ecting 
EMRIs with hitherto-unknown parameters (l24L sec- 
tion 3]. However, once EMRIs have been detected by 
more economical search algorithms f |l26l . section 3.1]; 
(194| ]]). precision modelling and matched filtering of the 
full LISA data stream become practical, allowing accu- 
rate measurements of the EMRI parameters, tests of gen- 
eral relativity, and other valuable astrop h ysica l measure- 
ments ( see, for example, jl25l Il95l4l99f : |l26l sections 4 
and 5]; [20q '). 

Gair [l27j | has recently updated past calculations |l24{ 
of LISA EMRI event rates and has calculated the red- 
shift z of the closest (zmin) and most distant (zmax) EM- 
RIs that LISA is likely to detect under a range of assump- 
tions about the LISA mission duration and hardware reli- 
ability, the supermassive black hole's spin, and t he EM RI 
rate per galaxy. For these calculations, Gair |127( as- 
sumed a detection threshold of /Otincsh — 30, where p is 
the EMRI signal-to-noise ratio after matched filtering. 
That is, Zmax is the redshift at which the strongest ex- 
pected LISA EMRI will have a signal-to-noise ratio (af- 
ter matched filtering) of /Othrcsh- Neglecting cosmological 
spacetimc curvature, the signal-to-noisc ratio for a given 
source scales inversely with z, so (neglecting Malmquist 
bias the signal-to- noise ratio of the closest LISA EMRI 
(which I take as an approximation to the strongest LISA 

EMRI) is thus Pmax ~ (2:max/2:min)Pthrosh- Table |T] givCS 

the resulting pmax for each of Gair's |l27l . table 4] LISA- 
performance and astrophysics assumptions. The Pmax 
values range from ~ 20 to as high as ~ 2000. 

In order to achieve these high signal-to-noise ratios, 
LISA will require matched filtering against accurate 
EMRI GW templates. In order to keep parameter- 
estimation errorfdj due to template inaccuracy below 
those due to statistical noise, the (template) EMRI GW 
phase must be modelled to an accuracy of A(p < C / pinax 
radians over the LISA mission lifespan, where the "de- 
generacy factor" C depends on the level of degeneracy 
between the different parameters for the particular anal- 
ysis being done. C is often estimated via the Fisher- 
matrix formalism (see, for example, [205l42ldj and refer- 



Malmquist bias is a selection effect in a brightness-limited sam- 
ple: nearby objects are included in the sample regardless of their 
intrinsic luminosity, but intrinsically-faint distant objects fall be- 
low the sample's minimum-brightness threshold and are thus 
omitted from the sample. The result is that the mea n intrinsi c 
luminosity of sample objects increases with distance |201| . l203l . 
In the present context, this results in the the 2:max EMRI be- 
ing intrinsically brighter than the z^\^ EMRI, which somewhat 
reduces pmax- 

These parameters might be those characterizing the EMRI sys- 
tem itself, those characterizing the deviation of the supermassive- 
body spacetimc from the Kerr metric, or those f or other t e sts o f 
general relativity (see, for example, [l25l Il95l , [l9^ . [l93, [l99| ; 
[m. sections 4 and 5]; [2o3, section 5]: [2041 ). 
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ences therein) 

LISA will observe an EMRI for TV ~ 27r • 10^ radians 
of GW phase (see, for example, [210l table I]), so the ac- 
curacy tolerance for the allowable GW phase error corre- 
sponds to a relative tolerance A(/)/(/) < C / {N p„^ax) for the 
instantaneous GW frequency. Table HIl gives these toler- 
ances for degeneracy parameters C = 1 (very optimistic) , 
C = 30 (reasonable for many tests-of-GR analyses), and 
C = 1000 (somewhat pessimistic). 

These GW error tolerances can be related to the re- 
quired accuracy in a self-forc e co mputation using the 
results of Huerta and Gair j210l table I], who esti- 
mate the effects of various 0{iJ,^) self- force effects - that 
is, 0{iJ? / ^ 10^^ fractional changes in the overall 
0(^) self force - on an EMRI's GW phase. They find 
that 0(^^) effects change the cumulative EMRI GW 
phase by '--^ 3 orbits (20 radians) over the ^ 10^-orbit 
LISA observation span. Equivalently, a 1 part per mil- 
lion (ppm) fractional change in the overall 0(^) self force 
changes the cumulative EMRI GW phase by approxi- 
mately 0.3 orbits (2 radians). The LISA EMRI phase 
error tolerances given in table |IT] thus correspond to the 
self- force accuracy tolerances given in table IIIll It's clear 
that self-force computations accurate to between roughly 
one part per million and one part per billion are required 



Signal-to-Noise Ratios of the Strongest LISA EMRIs 
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TABLE I. This table shows the estimated signal-to-noise ratio 
after matched filtering, pmax, of the closest (approximately 
the strongest) LISA EMRI sources, a is the dimensionless 
spin of the EMRI's central supermassive black hole, 7?.mw is 
the EMRI rate for the Milky Way galaxy, and "5yr,2d" and 
"2yr,ld" refer to different assumptions about the LISA mis- 
sion lifetime (5 versus 2 years) and hardware reliability (2chan 
= full configuration with 2 independent low-frequency inter- 
ferometer channels available; Ichan = degraded configuration 
with only 1 independent low-frequency interferometer chan- 
nel available). "*" marks values which are very uncertain due 
to small-A statistics in Gair's simulations [1271 . table 4]. 



Lindblom et al. [21ll421^il have carefully quantified a similar line 
of reasoning for the case of comparable-mass black hole binaries. 



to avoid degrading the parameter-estimation accuracy for 
the strongest LISA EMRIs. 

2. Self-Force Calculations 

As noted earlier, computing EMRI GW waveforms in a 
fully self-consistent manner requires calculating the met- 
ric perturbation induced by the particle - and the corre- 
sponding self- force - up to and including at l east 0(/x ^) 
terms f jl43l section 5.5.6], [lil section 11.1], fliliMg ]). 
but the theoretical formalism for doing this isn't fully de- 
veloped yet. 

However, in the near future some 0(/i^) effects are 
likely to be explored with "orbit correction" calcula- 
tions |l57L section 7]), where the 0(/ii) self force is used to 
calculate the time evolution of the orbit parameters. In 
order to reliably distinguish true 0(^^) effects due to the 
orbit correction from numerical errors in the 0(/i) self 
force, the 0(/x) self force needs to be calculated with a 
relative error <^ jj, ^ 10^^. 

This same argument should continue to hold once 
(if) future self-force calculations are able to include all 
0(/i^) effects and compute GW waveforms in a fully self- 
consistent manner. 

Ros enthal's work towards 0(/x^) self- force calcula- 
tions jl47l4l50l | suggests that the 0(/x) metric perturba- 
tion will be needed to high accuracy as an input into the 
0(/i^) calculations. 



Gravitational- Wave Phase 
Error Tolerance At/) (radians) 





C=l 


C = 30 C = 1000 


Pmax — 30 


0.03 


1 30 


Pmax = 300 


0.003 


0.1 3 


Pmax = 2000 


0.0005 


0.015 0.5 


Instantaneous 


Gravitational- Wave Frequency 


Fractional Error Tolerance /S.(f>/(f> 



C ^ 1 C = 30 C = 1000 

Pmax = 30 5 X 10^* 2 X 10-^' 5 X 10'^ 

Pmax = 300 5 X 10"^ 2 X 10^^ 5 X 10"'' 

Pmax = 2000 8 X 10"^" 2 X 10^* 8 x IQ-'^ 



TABLE II. This table shows the maximum errors allowed in 
an EMRI gravitational-wave template so that the resulting 
parameter-estimation errors for the strongest expected LISA 
EMRI do not exceed the statistical errors due to LISA's in- 
strumental and confusion noise levels, given various combina- 
tions of the EMRI signal-to-noise ratio pmax (after matched 
filtering) and the parameter degeneracy factor C. The error 
tolerances are expressed alternatively as a total phase error 
A(j> (radians), or as a (dimensionless) relative error in the 
instantaneous gravitational- wave frequency, A(j>/(j>- 
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Highly accurate self- force calculations are also valuable 
for helping to calibrate and constrain various terms in 
post-Newtoni an expansions multiple-body systems (see, 
for example, |l7ll . Il9ll l214j and references therein). 

Finally, highly accurate calculations of the 0(/i) self 
force are valuable as a test case for the intricate the- 
ory and computations involved. For the calculations re- 
ported here I use time-domain integrations of the metric- 
perturbation equations in the Barack-Ori mode-sum for- 
malism. In contrast, the most accurate published cal- 
culation of the self-force f or th is case, that of Detweiler, 
Messaritaki, and Whiting |l69l | , uses a frequency-domain 
approach with completely different numerical methods. 
Precisely because the two calculations are structured so 
differently, a verification of their agreement to high pre- 
cision serves as a useful check on both techniques and 
their respective theoretical formalisms 



Unfortunately, actual EMRI waveform calculations 
will likely be much slower than self-force calculations. 
For example, an orbit-correction calculation essentially 
requires time-integrating a set of coupled ODEs for 
the orbital-parameter evolution on radiation-reaction 
and longer timescales, with the ODEs' right-hand-side 
functions being given by a self-force computation |l57l . 
section 71) . Even the most efficient ODE-integration 
schemes |215| will require evaluating the right-hand-side 
functions (i.e., computing the self- force for some speci- 
fied intermediate orbit) hundreds of times in the course 
of a single orbit-correction calculation, so the need for the 
highest possible efficiency in the self-force computation is 
clear. 



III. SELF-FORCE CALCULATION VIA THE 
BARACK-ORI MODE-SUM REGULARIZATION 



B. The Importance of High Efficiency 

The precision modelling and matched filtering of a 
single already-detected EMRI is essentially a many- 
parameter nonlinear least-squares fitting process, and 
thus requires generating many trial waveforms. More- 
over, this process should be repeated for each stro ng 
EMRI source, of which there will likely be many 
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With current methods, a single EMRI self- force calcu- 
lation takes between one -half and one cpu-week at the 
10""* relative-error level jl77l section III.E]. This is al- 
ready unpleasantly slow, and raising the accuracy to the 
< 10~^ relative-error level will slow the computation by 
another factor of ~ 100 although parallelization should 
be easy. 

Self-Force Relative Error Tolerance 





C = 1 


C = 30 


C = 1000 


pmax 


= 30 2 X 10"* 


5 X 10"^ 


2 X 10"^ 


pmax 


= 300 2 X 10"^ 


5 X 10"** 


2 X 10"** 


pmax 


= 2000 3 X 10"^° 


8 X 10"^ 


3 X lO"'' 



TABLE III. This table shows the maximum relative errors al- 
lowed in an EMRI self-force computation so that the resulting 
parameter-estimation errors for the strongest expected LISA 
EMRI do not exceed the statistical errors due to LISA's in- 
strumental and confusion noise levels, given diflFerent combi- 
nations of the EMRI signal-to-noise ratio pmax (after matched 
filtering) and the parameter degeneracy factor C. 



Sago, Barack, and Dctwcilcr Il62l and Barack and Sago Il77ll 
have previously compared time- and frequency-domain self-force 
calculations. Comparisons of self-force calculations w ith post- 
Newtonian expansions (see, for example, llTlillQll , |214|| and ref- 
erences therein) also implicitly ch eck the correctness of both. 
Like my code, Barack and Sago's IITTH code uses globally 4th or- 
der finite difi^erencing, so a X 100 accuracy improvement requires 



In this section I briefly outline the Barack-Ori mode- 
sum regularization procedure for computing the self 
force, for the special case of a scalar particle in a cir- 
cular geodesic orbit in Schwarzschild spacetime. A more 
detailed account can be found in the original works by 
Barack and Ori jl63l - ll67j . I defer most discussion of nu- 
merical methods for this calculation to section ITVl 



A. Schwarzschild spacetime 

Consider Schwarzschild spacetime of mass M, and in- 
troduce ingoing and outgoing null coordinates u and v 
respectively, so the line element is 



ds^ = -f{r) dudv + r^{dO^ + 



sm 



(1) 



where r is the usual areal radial coordinate, /(r) = 1 — 
2M/r, and {B^Lp) are the usual polar spherical angular 
coordinates on a 2-sphere of constant r. It's also useful to 
define the Schwarzschild time coordinate tschw = \{v+u) 
and the "tortise" radial coordinate 



2Mlog 



r 
2M 



1 



(2) 



It's convenient to define the specific energy £, specific 
angular momentum L, and orbital frequency a; of a test 
particle in a circular geodesic orbit at the areal radius r, 

f{r) 



£{r) 
L(r) 




- 3M/r 



(3) 
(4) 
(5) 



a X \/lO increase in resolution, which costs a factor of 10 in CPU 
time for a l-|-l-dimensional evolution. 
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B. The Scalar Field 



I take the real scalar field $ to satisfy the equation 



■dr, 



(6) 



where q is the particle's scalar charge and r is proper time 
along the particle's worldline Xp = Xp{T). Specializing to 
the particle being in a circular geodesic orbit at areal 
radius r = r^, aligning the equator of the coordinate 
system (0 = |) with the particle orbit, and changing the 
variable of integration from proper time r to coordinate 
time ischw, this becomes 



□$ = -1^^5(r - rp)5{e - f )<5(^ - c^ptgchw) , (7) 

fp Op 



where (and henceforth) the subscript "p" denotes evalu- 
ation on the particle's worldline r = rp. 



Now expand r$ in spherical harmonics {YimiO,^)} 
(with normalization given by (|12|) below) by defining the 
complex scalar fields = (l>eni{tschw, r) such that 

oo e 

e=0 rn=~l 

(8) 

Each satisfies the inhomogeneous linear wave 

equation 

□0Cm + (r)0£„ = 5'fm(tsciiw)'5(r - Vp) , (9) 
where the potential Vi and source term Sim are given by 



Schw) 



4 



2M ^ £{£+!) 



(10) 



TpSp 



■ exp(-ima;ptsciiw) , (H) 



with the (real) coefficients {aim} defined by 



y,„(0=f,(p) =a,„e™'^, (12a) 



I.e. 



, ,,(i+m)/2 /2£+l /(£ + m-l)!!(£-m-l)!! ^ 

{-Ir \ \ m rrrh ttt-^ if ^-m is even f.^,. 

^ ^ V 47r Y {i + m)\l {£ - my.l , (12b) 

if £—m is odd 



where the "double factorial" function is defined by 



n-(7i-2)!! ifn>2 
1 if n < 1 



(12c) 



Each (j)im can be obtained by numerically solving the wave equation ([9]). I discuss the problem domain and boundary 
conditions for this equation in section fill D[ and I discuss the numerical solution in section flV Al 



C. Computing the Self-Force 

Assuming that the complex scalar field (f>im is known for each (^, m) , the contravaraint radial component -Fgoif of 
the 0(/i) self force may be computed as described by Barack and Sago [189|: For each i > 0, define 



d{(j)im/r) 



dr 



(13) 



tschw,''=r* 



where r = refers to computing the one-sided derivative as r approaches the particle worldline either from the 
outside (-I-) or the inside (— ), in both cases on a slice of constant tschw For finite-differencing purposes, it's convenient 
to transform this derivative into one with respect to r*: since dr^jdr = 1/ f , we have that 

d{(i)lm/r) _ 1 d(j)lm _ (km ^^^^ 

dr fr dr^ 

Now (following Barack and Lousto [l73l |) observe that under the transformation m — — m, the wave equation's 
source term Sim defined by (|lip transforms to its complex conjugate. Since the wave equation's potential Vi is real 
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and independent of m, this means that the equation's solution (jjgrn also transforms to its complex conjugate. Thus 
(using (|12ap ). simplifies to 



Ff^ (ischw) = ^ aem exp{imu}ptschw) 



d{(j)im/r) 



m=0 

where for any quantities X^rm we define the notation 



dr 



(15a) 



tSchw,''=''p 



(15b) 



m=0 



Following Barack and Ori [166| . the contravariant ra- 
dial component of the self-force at any point on the par- 
ticle's worldline is then given by 



Schw ) 7 



(16a) 



£=0 



where the regularized self-force modes F^^^^ are given by 

i^Sgfechw) = Hischw)T(^+5)A(rp)-B(rp) , (16b) 

where (for a particle in a circular geodesic orbit in 
Schwarzschild spacetime) the regularization coefficients 
A{r) and B{r) are given by 



Mr) = 4^ 



B{r) 



q^£^[E{w)-2K{w)] 



.2 7r/V3/2 
where V and w are given by 

V = 1 +LVr^ 



(17a) 
(17b) 

(18) 
(19) 



and K{w) and E{w) are the complete elliptic integrals of 
the first and second kinds respectively, 

K{w) = / ^ da; (20a) 



1 

\/ 1 — w sin" X 

/•V2 

-E(zi;) = / V 1 — w sin^ x dx . 
Jo 



(20b) 



Barack and Ori |l66l | have shown that F^^^^ ~ F^ 

and hence that F^^f = F^^J . In view of this the 
and superscripts may be dropped, and we may 
rewrite (jl6ap as 



^self = ^ ^£,re 



(21) 



without ambiguity. However, for numerical purposes 



it's still very useful to compute both expressions F^ 
~) . 

rog' 



and F^ 2 ■ I discuss this in section HV El 



D. Problem Domain and Boundary Conditions 

The wave equation (|9]) is naturally posed on an in- 
finitely large domain with boundary conditions at infin- 
ity appropriate for an isolated system in an asymptot- 
ically flat spacetime. However, for numerical purposes 
it's convenient to inst ead f ollow an approach suggested 
by Barack and Lousto |173| , solving ® on a large but fi- 
nite domain using arbitrary initial data and/or boundary 
conditions. These introduce a burst of spurious "radia- 
tion" dynamics into the solution , but fortunately this 
spurious radiation dies out quite quickly as one moves 
away from the initial slice(s) and/or the problem-domain 
boundariesEl The self-force is defined along the parti- 
cle's worldline, and its value at a given event Q on that 
worldline depends only on and V4'im at Q. The 
effect of the spurious radiation can thus be made negli- 
gible by choosing a sufficiently large numerical problem 
domain whose initial slice and/or boundaries are suffi- 
ciently distant from Q. 



E. The Tail Sum 

The definition (PT|) of the self-force F^*^ is written in 
terms of an infinite sum X)£=o regularized self-force 
modes i^f,rog- For numerical purposes a fin ite expression 
is needed. Following Barack and Sago |l89l . section III.E] , 
partition the infinite sum (|2ip into a finite "numerical 
force" sum of the modes with i < K and an infinite "tail 
force" sum of the modes with £ > K' = K+1, where 
if ~ 30 is a numerical parameter: 

-Ffsolt = i^solf.num + i^solf.tail (22a) 
K 

(22b) 



oo 



seli,num — / , ^ £,rcg 
oo 

^sclf,tail 



(22c) 



I have seen no ev idence of the Jost "persistent junk" solutions 
discussed by |216|I . 
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Once the regularized self-force modes Fi^, 



for < ^ < if, the numerical force K 



scll,num 



are known 
is easy to 

compute from the definition (j22bp . The tail force i^soif,taii 
can be es timated using the known large-£ series expan- 
sion jl69l . equation (12)] 



p even 



(23) 



where the {cp} are coefficients not depending on and 
the basis functions fp{i) = 0{£~p) are given by 



f2i£) 

m 



1 

[i - m - m - m - m + m + m + m + d 



(24a) 
(24b) 
(24c) 
(24d) 



Typically only a few terms in this series are needed to give an excellent approximation to F^,reg- 

For a particle in a circular geodesic orbit in Schwarzschild spacetime, Detweiler, Messaritaki, and Whiting |16 
have shown that the coefficient C2 is given by 



C2 




2r2(rp-2Af) 



3M 



M{rp - 2M) ^r,-M){rp-iM) ^ 

aA/^ _oA.r^ "^1/2 



3r4(rp - 2M) 



(rp - 3Af)(5r2 - IvpM - UA'P) 
16r4(rp-2A/)2 



G 



3/2 



3{rp-3M)^{rp + M) 
16r4(rp -2M)2 

I 



G 



5/2 



(25) 



where the leading factor of —1/4 converts from the nor- 
malization used by Detweiler, Messaritaki, and Whiting 
to that used here, and where Gp (a special case of a Gauss 
hypergeometric function) is given by 



2 f^/'^ 

Gp^- (l-asin^a;)-^ 
I" Jo 



able subset of the numerically-computed -Fi,reg values. I 
discuss the numerical computation of this "tail fit" in 
section IIVDI 

Once the {cp} coefficients are known, the tail 
force (|22cp is then given by 



dx , 



(26) 



sclf.tail 



with a = M/{rp — 2M). G±i/2 can also be written in 
terms of the complete elliptic integrals ([20)) . 



oo 

E 

e=K' 



p even 
p>2 



(28) 



where 



G_i/2 = -E{a) 
Gi/2 = -Kia) . 



(27a) 
(27b) 



(29) 



=K' 



The C4 and higher coefficients aren't known analyti- 
cally, but they can be estimated numerically by least- 
squares fitting the tail-series expansion (j23|) to some suit- 



Using the Maple symbolic algebra system f j217l |. 
|http : //www .map lesof t . com/) version 11) to evaluate 
the sums |29|) F^ I find that the first few Tp are given 

by 



^® These sums can also be evaluated by hand by first using partial 
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r2 
r4 



K' 

{K'-\){K' + \) 

K[ 

3(A''-f)(/v'-i)(X' + i)(if' + f) 

IC_ 

5(if' - |)(A- - ' - i)(A- + \){K' + + I) 



7{K' - |)(A- - - - + UK' + UK' + f )(/^' + I) 



(30a) 
(30b) 
(30c) 
(30d) 



IV. NUMERICAL COMPUTATION OF THE 
SELF-FORCE 

In this section I describe the numerical methods I use 
for high-accuracy self-force calculations. 



A. Numerical Solution of the Wave Equation dS]) 

1. General Numerical Scheme 

Near the particle worldline the complex scalar field 
(j)irn has 0(1) amplitude and rapidly oscillating phase in 
both space and time, but the field amplitude decreases 
quickly with increasing distance from the particle world- 
line. This high dynamic range suggests the use of a 
mesh-refinement method to resolve the fast oscillations 
without the computational cost of of maintaining this 
high resolution everywhere in the numerical domain. The 
numerical method also needs to accommodate the non- 
differentiability of (pirn across the particle worldline. 

To avoid the numerical complications of ex plicit 
boundary conditions, I follow Barack and Lousto jl73j 
and use a characteristic (double-null) numerical evolu- 
tion scheme, with a "diamond-shaped" problem domain 
which is a square in the characteristic variables u and 
■y, {u,v) G [wi„in,Mmax] X [umin,Wmax]. With this domain 
the (arbitrary) initial data (pim = is applied on the 
"southwest" and "southeast" grid faces v = Wmin and 
u = Mmin respectively; I place the domain such that the 
particle worldline r = Vp symmetrically bisects the do- 
main. Figure [T] illustrates the problem domain and par- 
ticle worldline. This type of problem setup has been used 
successfully for a number of othe r sel f -forc e calculations, 
including (for example) those of jl76l Il77( . 

To numerically solve the wave equation © on this do- 
main, I use a characteristic adaptive mesh refinement 
(AMR) numerical scheme with 4th order global finite dif- 
ferencing accuracy. I have described this scheme in detail 



elsewhere |l78l |. Briefiy, the underlying (unigrid) finite 
differencing is a standard double-null di amond-in tegral 
scheme with square grid cells in {v,u) ex- 
tended to provide globally 4th order finite differencing 
accuracy in a manner similar to that of [l^76, 223J. The 
AMR algorithm it is very simila r to the stand ard Cauch y 
Berger-Oliger AMR algorithm ([Ull; see also j226l4230j ). 
slightly modified as suggested by Hamade and Stew- 
art [23ll | to accommodate the characteristic evolution. 
The AMR algorithm treats treats w as a "time" coordi- 
nate and u as a "space" coordinate: the evolution inte- 
grates V = constant slices successively in the direction 
of increasing v, with each slice completely integrated (in 
the direction of increasing u) before the integration of the 
next slice begins. 

The AMR algorithm begins with a relatively coarse 
"base" grid which covers the entire problem domain; dur- 
ing the evolution the algorithm dynamically (adaptively) 
constructs a hierarchy of finer "child" grids, each a factor 



'Schw 




fractions, after which each sum telescopes, then finally undoing 
the partial fractions to further simplify the result. I have explic- 
itly verified ll30at . llSObt . and l30ct in this way. 



FIG. 1. This figure shows the overall problem domain, and 
the {u,v) and (ischw,*"*) coordinates. The vertical line marks 
the particle worldline. Mesh refinement is inhibited in the 
V-shaped shaded region lOOM wide bordering the "south- 
east" and "southwest" grid faces. The self-force is measured 
along the region of the particle worldline marked by horizontal 
hatching. 
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of 2 finer than, and spatially nested inside, its "parent" 
grid. The fine grids typically only cover small subsets of 
the problem domain. 

The AMR algorithm is controlled (the " adaptive" part 
of AMR) by comparing an estimate A of the numerical 
solution's local truncation error (LTE|13 with aspeci- 
fied tolerance eito. If (after median smoothing (l78l sec- 
tion 4.2]) A > eitc then the algorithm adds another level 
of mesh-refinement to better resolve the solution. 

As well as AMR, the numerical scheme also uses fixed 
mesh refinement (FMR): following [230t . my AMR code 
has options to record the placement and grid spacing of 
each refinement level generated by the AMR algorithm. 
This can then be "played back" with each grid refined by 
a chosen small- integer factor iVfmr. FMR is useful both 
for convergence tests, and in some cases for circumventing 
floating-point roundoff limits on my AMR scheme (these 
are discussed in section lTV A2[) . 

Because of the characteristic evolution scheme, the lo- 
cal finite difi^erencing must actually be 6th order acc urate 
in order to achieve a 4th order global accuracy (see [itI 
section 3.1] and references therein). Similarly, the global 
accuracy generally scales as s^l^ , or e^l^^^ if FMR is used, 
where the "effective error tolerance" is eitc,ofi' = ^^tc/Nf^^ 
(cf. discussion in section IVDl particular figure [TT|). 

2. Extended Floating-Point Precision 

Floating-point numbers are only represented and com- 
puted with finite accuracy; typically each floating-point 
operation introduces a small roundoff error of fractional 
size <efp, where Efp = 2'^'^ « 2.2x10-^^ for IEEE- 
standard double-precision floating-point arithmetic 

There are (at least) two different parts of my numerical 
scheme for solving the wave equation Q which may be 
limited in accuracy by floating-point roundoff effects: 

|1] The first and most obvious way in which floating- 
point roundoff effects limits the achievable accu- 
racy of the numerical scheme is the finite-difference 
computation of (^fm at each successive gr id p oint. 
This computation (described in detail in |l78l ap- 
pendix A. 2]) involves ~ 50 floating-point opera- 
tions. In the absence of fortuituous error cancel- 
lations, this computation contributes a relative er- 
ror of (TEfp at each grid point, where cr > 1 reflects 



The LTE is a measure of the local accuracy with which the finite 
difference equations approximate the underlying PDE (here the 
wave equation ) . More precisely, the LTE is a pointwise norm 
of the discrepancy that would result if the exact solution of the 
PDE were substitu ted into the finite difference equations at a 
grid point [23^^231 . 

More precisely, £fp, usually known as the "machine epsilon", is 
defined as the smallest positive floating-point number such that 
1 © £fp 7^ Ij where © is the floating-poi nt ad dition oper ator . 
This is discusse d in d etail by, for example, |236| . chapter 2]; |237| , 
chapter 2]; and [238ll and references therein. 



the error-propagation properties of the computa- 
tion (which I have not analyzed in detail). 

|2] The second way in which floating-point roundoff 
may limit the achievable precision of my numer- 
ical scheme is via the AMR algorithm: My code 
estimates the LTE by comparing the standard nu- 
merical computation of (/> at a grid point with an al- 
ternate lower-resolution computation which spans 
the most recent 2 grid points in v and u with a 
single finite differencing step |l78l . equation (6)]. 
If the difference between computed in these two 
ways isn't well resolved by the floating-point arith- 
metic, the LTE estimate will be unreliablelil In 
practice, taking into account the normalization fac- 
tors in the actual LTE estimate, I ensure reliable 
operation of the AMR algorithm by limiting it to 
an LTE-cstimate tolerance ^ Sfp. 

One way to circu mvent the AMR LTE-estimate 
limit |2| is (following j230l |) to record the placement 
and grid spacing of each reflnement level generated by 
the AMR algorithm, then "play back" this with each 
grid refined by a chosen small-integer factor Afmr- This 
"flxed mesh refinement" (FMR) reduces the global finite- 
difference truncation error (the cumulative effects of the 
LTE in all the grid cells in the entire nu meric al integra- 
tion) by very close to a factor of Nf^^ |l78l . figure 7], 
at the cost of an increase in the code's running time by 
a factor of Nf^^^.. However, the per-grid point rounding 
error limit |1] cannot be circumvented in this way. (In 
fact, FMR may worsen the overall floating-point round- 
off errors in the self-force by a factor of Nf^-^^ or more due 
to the larger number of individually-smaller grid cells in 
the integration.) 

Due to the AMR LTE-estimate limit |2] , I restrict the 
AMR algorithm to a tolerance £\tc ^ 10~^^ when us- 
ing standard IEEE double-precision floating-point arith- 
metic. FMR can improve this considerably, but beyond 
A'^finr ~ 6, the per-grid-point rounding error limit [IJ be- 
comes increasingly severe, and my error estimates for the 
individual -Fi,regi the tail flt, and the overall self-force all 
become less reliable. 

To further investigate the effects of floating-point 
rounding errors in the numerical solution of the wave 
equation I extended-precision floating-point arith- 
metic. In particular, on Intel x86 and compatible proces- 
sors my AMR code for solving the wave equation ^ 
can optionally use IEEE "double-extended" floating- 
point arithmetic (typically specified in C/C++ as "long 



If A is unreliable, then (even after the smoothing) we might well 
have A > en^ somewhere on each new slice, no matter how small 
the grid spacing. This would cause the AMR, algorithm to effec- 
tively infinite-loop, continually adding further reflnement levels 
until it runs out of memory. Although a limit on the maximum 
refinement level could prevent this, the algorithm would still be 
refining inappropriately, causing the computation to be very in- 
efficient. 
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double"). This provides a relative accuracy of £[p = 
2-^3 I I X 10-19, a factor of 2" = 2048 times more 
accurate than IEEE double precision. This lowers the 
AMR LTE-estimate hmit {2] to enc > 10^^^, with only a 
modest performance penalty compared to standard IEEE 
double precision (at the same accuracy setting my code is 
about a factor of 2 slower in long-double than in double 
precision). 

Note that even when using extended-precision 
arithmetic in this way, once the gradients 
d{(j)em/r) /dr\ _ ± wee known along the particle 
worldlinc, the remainder of the self-force computation 
is considerably less sensitive to the floating-point 
arithmetic precision. I thus use standard IEEE double 
precision for computing each regularized self- force Fp^^^^ 
the numerical force (j22bp . the tail fit and tail force 
and the error estimates for these quantities!^ 



B. Parallel Execution 

Even with AMR, self-force computations are still very 
expensive, so it's useful to parallelize them as much as 
possible. Fortunately, the self-force problem is trivially 
parallelizable by distributing the solution of the wave 
equation ((9)) to different processors for different (£, m). 
Because no communication is needed between the com- 
putations for different {^, m), this requires very little com- 
munications bandwidth, and overall performance scales 
almost linearly with the number of processors usedF^ 

For the results presented here, I used between 10 and 
15 processors of a local workstation cluster, with a shared 
NFS file system to collect the results from each proces- 
sor's computations. Each processor was either a 2.5 GHz, 
2.8 GHz, or 3.2 GHz Pentium 4. 



C. Regularization Coefficients 

I compute the regularization coefficients A{r) and 
B{r) from the definitions ((T7)) . evaluating the K 
and E complete elliptic integrals ([^0)1 using the 
e llpk and ellpe subroutines from the Cephes library 
([240[, http : / /w ww . net lib . org/ cephes, release 2.2 
dated July 1992). 



To (slightly) red uce fl o ating -point roundofi^ errors, I use Ka- 
han summation ( |239|I : l238l theorem 8]) when evaluating the 
sums ifTsjl . I|16a|l and lf28)l . 

In the parallel-computing community, this type of problem is 
known as "embarrassingly parallel" , in the sense that it's such 
an easy test case for parallel hardware that one should be em- 
barrassed to report parallel-speedup results for it. 



D. The Tail Fit 

I consider two cases for the tail fit: 

• For the most accurate computation possible (as- 
suming the particle to be in a circular geodesic 
orbit), I compute the tail- fit coefficient C2 from 
the expression ([25]) . evaluating each Gp via di- 
rect numerical integration of the definition ([26|), 
using the dqags subroutine (revision date 198 3 
May 18) from the QuADPACK library f j241 
http: //www. net lib. or g/quadpack) 



As noted in section UlI E[ the C4 and higher tail- fit 
coefficients can be estimated numerically by least- 
squares fitting the series expansion ((23|) to some 
suitable subset of the numerically-computed -F^,reg 
values. For the accuracies obtained in this paper, 
it suffices to keep only terms up to and including 
the 0(£^^) term, so the tail fit only includes the 
coefficients {04, cg}. 

• Alternatively, to simulate the accuracy to be ex- 
pected for a particle in a generic non-circular orbit 
(where the C2 coefficient isn't known analytically 
for the form of the mode-sum regularization used 
here) 1 also consider the case where C2 is included 
in the tail fit, i.e., where the coefficients {c2, C4, ce} 
are fitted simultaneously. 

Whichever set of coefficients are fitted, computing the 
tail fit numerically requires some care, because the ba- 
sis functions {/p} defined by p4)) are nearly degenerate 
(linearly dependent), causing the tail fit to be quite ill- 
conditioned. That is, there are linear combinations of 
the basis functions bpfp where the linear-combination 
coefficients {hp] have unit 2-norm (call these "unit- 
coefficient- norm" linear combinations), yet the linear 
combination bpfp is very small relative to the largest 
of the {/p}. The fitted coefficients {cp} are relatively 
uncertain in the direction of any such {6p}, which intro- 
duces additional uncertainty into the tail force i^soif,taii 
computed via ([28]) . 

Figure [2] illustrates the near-degeneracy of the {/p}, 
showing very small unit-coefRcient-norm linear combina- 
tions of various subsets of the {/p}, and table IIVI gives 
the corresponding condition numbers The condi- 

tion numbers are primarily determined by how many 



I have also explicitly verified that the identities (I27t hold to very 
high accuracy (a few parts in 10^^) for my numerical implemen- 
tation. 

Haas and Poisson jl70| | have computed the equivalent of the 
C2 coe fficie nt for a different form of mode-sum regularization, and 
Haas |176|| has used this successfully in a numerical computation 
of the self-force on a scalar particle in a generic (non-circular) 
orbit in Schwarzschild spacetime. 
^■^ These very-small linear combinations can be determined from 
a singular value dec o mpos ition (SVD) of th e least-square s fit' s 
design matrix ( |242|| ; [23^ , chapter 9]; [2371 . section 6.8]; [243| : 
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coefficients {ck} are fit simultaneously: fitting 2, 3, or 
4 coefficients simultaneously gives a condition number of 
of K ^ lO'', 10^, or 10^ respectively. 

Because of this ill-conditioning, it's much better (yields 
more accurate results) to perform the tail fit using a QR 
or singular valu e dec omposition, rather than via the nor- 
mal equations ( j245l section 2.2]). (If the normal equa- 
tions were used, the effective condition number would be 
roughly the square of that given here, thus greatly in- 
creasing the effects of floating-point roundoff errors on 
the results.) 

Much of this ill-conditioning is due to the widely dif- 
fering magnitudes of the different basis functions (this 
can be seen in figure [2|) , and can be greatly alleviated by 
simply renormalizing the basis functions to have similar 
magnitudes over the range of i used in the tail fit. To 
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this end, I define 




(31) 



where the parameter £ (taken here to be 20) is the £ at 
which all the normalized basis functions will now have 
unit magnitude. Table ITVl also gives the condition num- 
ber for fits using various subsets of the normalized basis 
functions {/p}. The normalized basis sets have much 
smaller condition numbers, and correspondingly lead to 
significantly smaller floating-point roundoff effects in the 
tail fits. 

I compute the fitted coefficients {cp} and their covari- 
ance matrix using the p;sl_multif it_wlinear _svd sub- 
routine from the GNU Scientific Library f [246l |. ver- 
sion 1.12), using the normalized basis functions {fp{£)} 
defined by ([3T|l F^ I assign the individual data points the 
weights {SF^:"Y'. 




ell 

— f2(ell) • linear combination of t^Jg 

— f4(ell) + linear combination of i2'U'^6 

— iQ{e\\) X linear combination of f2,f4,f6,f8 
f8(eii) 

FIG. 2. [Color online] This figure shows the tail-fit basis 
functions /2, fi, fa, and fs, together with 3 very small unit- 
coefficient-norm linear combinations of the basis functions. 



[244l . section 15.4]; [245j V For present purposes the condition 
number k can be interpreted as the ratio of the largest 2-norm of 
any basis function to the smallest 2-norm of any unit-coefficient- 
norm linear combination bpfp of the basis functions. Thus 
1 < K < oo, with K = 1 describing an orthonormal basis set, 
K 2> 1 describing a nearly degenerate basis set, and k = oo 
describing a perfectly degenerate (linearly dependent) basis set. 
Small errors in the input data and/or computation of the fit - 
including in particular floating-point roundoff errors — are am- 
plified by a factor proportional to k in the outputs of the fit (the 
fitted coefficients {cp}), and thus also in the tail force Fg^if tail 
computed via I I28I1 . 



E. Internal Error Estimates 

I now consider the numerical computation of error es- 
timates (bounds) for the individual regularized self-force 
modes, the numerical force, the tail force, and the overall 
self-force. In this section I consider only "internal" error 
estimates, those which can be computed from (as part of) 
a single self- force calculation. In section HV Fl I consider 



Coefficients 


Condition Number k 


Being Fitted 


Basis is {fp} 


Basis is {fj 


{c4, ce} 


2.6 X 10^ 


11 


{C4, C6, Cg} 


5.9 X 10^ 


100 


{C2, C4} 


2.2 X 10^ 


8.3 


{c2, C4, Ce} 


4.4 X 10*^ 


66 


{C2, C4, Cg, Cg} 


8.0 X 10^ 


460 



TABLE IV. This table shows the condition number n of the 
tail fit (more precisely, of the fit's design matrix if all the 
regularized self-force modes -Fi.reg are taken to have unit un- 
certanties) when various sets of tail-fit coefficients {cp} are 
fitted and either the {fp} or {fp} basis functions are used in 
the fit. The Fe. rog are assumed to be given for 20 < ^ < 30, 
£ = 35, and £ = 40, as is the case for the numerical results 
presented in section [V] However, the condition numbers de- 
pend only relatively weakly on the precise set of £ used and 
the relative uncertainties of the different Fi^j-^g. 



gsljnultlf it_wlinear_svd actually does its own scaling inter- 
nally, similar to my normalization I I31I1 . However, not all QR- 
or SVD-based least-squares fitting routines do this; for example, 
the widely-used SVD-based routines given by |244| . section 15.4] 
do not perform such an scaling internally. 
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"record-playback" error estimates derived from compar- 
isons between a low- and high-accuracy pair of self-force 
calculations, and in section HV Gl I consider "actual" er- 
rors derived from comparisons between a self-force cal- 
culation and a different and much more accurate calcu- 
lation. 



Individual Regularized Self-Force Modes Fi^r 



As noted in section llil CI F^^l^ — -F/7og- However, due 
to the finite-difference truncation errors in the numeri- 
cal solution of the wave equation ([3]), the numerically- 
computed values of Ff~ll^ and Ff~}„ will differ slightly. I 
use this difference to derive an error estimate (more ac- 
curately, an error hound) for each regularized self-force 
mode, 



— 2 \fi,r 
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F, 



{-) 
'.rcg 
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(32) 



[Notice that this internal error estimate does not de- 
pend in any way on the use of an AMR algorithm to solve 

the wave equation ([9]): i^/^jg ^^"^ -^i^ig would both still 
be well-defined even in a unigrid simulation, and their dif- 
ference would still be a measure of the finite-differencing 
errors. 

In section fV C II I present numerical evidence that these 
error estimates provide reasonable (in fact, somewhat 
conservative) estimates of the actual numerical errors in 
the individual regularized self-force modes. 



2. The Numerical Force Fseif,num. 

The propagation of the individual regularized self- force 
modes' error estimates (|32|) through the numerical-force 
computation (j22b[) is non-trivial, because we don't a pri- 
ori know whether or to what extent the actual errors 
in the individual regularized self-force modes F^.rog for 
different ^ are correlated. (Since all the modes are calcu- 
lated using the same basic numerical scheme, some degree 
of correlation in their errors would not be implausible.) 

To investigate this question, I consider two extreme 
cases for how the numerical force's error estimate 



5F^ 



(internal) 

3lf 



might be defined: 



• If the actual numerical errors in different modes 
are statistically independent, then the individual 
modes' error estimates should be added in 

quadrature. 
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(33a) 



error estimates ((32)) should be added arithmetically. 
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6F^ 



(internal) 
;lf,num 



(internal) 
eg 



(33b) 



i=0 



Barack and Sago |177l | use this formula to com- 
pute the numerical-force error given the (record- 
playback) error estimates of the individual modes. 



In section rVC 21 1 present numerical evidence that the 
actual numerical errors in different modes are, if not com- 
pletely independent, then at least weakly enough corre- 
lated that the quadrature sum (|33ap provides a reliable 
error estimate for the numerical force, whereas the arith- 
metic sum (j33bp systematically overestimates the errors 
in the numerical force by a factor of ~ 3. 



3. The Tail Force Fseif,taii 

For fitting the large-^ series expansion and com- 
puting the self- force tail forceFscif,taii via p8)) . the same 
issue of statistical independence versus correlation of the 
actual errors in i^£,reg for different i arises again: 

• If the actual numerical errors in different modes 
are statistically independent (and the individual 
modes' error estimates ([32]) are treated as the stan- 
dard deviations of Gaussian distributions) , then the 
standard the ory o f liii ea r lea st-squares fitting can 
be applied {^M, HH; HI section 15.4]). The 
tail fit then provides the covariance matrix C for 
the fitted coefficients {cp}l^ Since the tail force 
can be written as the linear combination ([25]) of 
the fitted coefficients {cp} with linear-combination 
coefficients {Fp}, it's then easy to compute the es- 
timated uncertainty in the tail force. 



'SF 



(internal) 



self.tail 



EC: 



F F 



(34a) 



where the sum is over only the linear-combination 
coefficients {Fp} corresponding to the fitted coeffi- 
cients {Cp}. 

• Alternatively, I can approximate the worst-case er- 
rors in the tail force (|28p without assuming any- 
thing about the statistical independence of the ac- 
tual numerical errors in different modes, as fol- 
lows: Suppose the set of £ for which Fg^^eg is 
used in the tail fit is ^2, ^3, ■ • ■ , ^q}- For each 
fc g {1, 2, 3, . . . , Q}, suppose r]k G {—1, 0, -1-1}, and 
define F^ = F,,,^, + ry. JF^^^^^^^ Then for 
each of the 3'^ possible combinations of 771, 7^2, 



• If the actual numerical errors in different modes 
are perfectly correlated, then the individual modes' 



If the statistical assumptions don't actually hold, C is more ac- 
curately termed the formal covariance matrix. 
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r^Q, I perform a separate "trail" tail fit of the self- force mode -F^^rcg by 



-,(trial) 



and com- 



series expansion (f23|) to all the 

pute the corresponding tail forceFjg["taii 
Finally, I take the extreme range of these tail forces 
^scif'taii among all 3*^ trial tail fits as a worst-case 
error estimate (^i^J^JJ^'^^j]^'^ for the tail force Fseif,taii, 



(internal) 

;lf,tail = 



max 



p(trial) _ p 
^self,tail -f^sclf.tail 



(34b) 



In section rVC 31 1 present numerical evidence that the 
statistical error estimate p4ap is moderately conserva- 
tive, overestimating the actual numerical errors in the 
tail force i^scif.taii by a factor of while the worst-case 
error estimate p4bp overestimates the actual errors by a 
factor of 5. 



4- Overall Self-Force 



Given the internal error estimates p3ap and (|34ap for 
the numerical force and tail force respectively, the ques- 
tion of their statistical independence or lack thereof arises 
once again when computing the overall self-force itself 
via (|22ap . I again consider two possible choices for an 

internal error estimate SF^ 



(internal) 
self 



for K 



self- 



• As a best-case estimate (errors perfectly indepen- 
dent), I take the quadrature sum 



SR 



(internal) 



sclf.num 



-5^1;.;:^) (35a) 



• As a worst-case estimate (errors perfectly corre- 
lated), I take the arithmetic sum 



(internal) r 77!(internal) 

self,num 



self 



r 77!(internal) 
" "^self.tail 



(35b) 



In section IV C 41 I present numerical evidence that 
while both of these error estimates are fairly reliable, 
the arithmetic-sum error estimate (j35b[) tends to give a 
slightly more accurate estimate of the actual errors than 
the quadrature-sum error estimate pSap . 



F. Record-Playback Error Estimates 

To validate the internal error estimates, I pair each 
"record" AMR solution of the wave equation Q with a 
corresponding "playback2" numerical solution incorpo- 
rating FMR of the recorded grid structure by a factor 
of Nf^r = 2, in the manner discussed in section IIVA2I 
My numerical code shows excellent 4th order conver- 
gence |l78l figure 7], so the finite-difference truncation 
errors in each playback2 evolution are very close to a fac- 
tor of -/Vf^jr = 16 smaller than those of the correspond- 
ing AMR record evolution. I thus define the "record- 
playback" error estimate for each "record" regularized 
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r^(record) 
" ^,reg 



-i(playback2) 



(36) 



[Here I'm implicitly assuming that finite-difference trun- 
cation errors are the major contributor to the overall er- 
ror in each F£,reg- As discussed in section IIV A21 this is 
true in practice so long as the AMR error tolerance eitc 
isn't too small, cf. discussion in section IV Dl ] 

This same record-playback technique can also be used 
to estimate the numerical errors in the "record" nu- 
merical force i^soif.num, tail force i^scif.taii, and total self- 
force itself, 
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(37a) 
(37b) 
(37c) 



Notice that unlike the internal error estimates dis- 
cussed in section IIV El which can be computed from a 
single numerical solution of the wave equation the 
computation of any of these record-playback error esti- 
mates requires a pair of numerical solutions of different 
accuracy (in this case, record and playback2); the record- 
playback error estimate is only computed for the lower- 
accuracy (in this case, record) member of the pair. 



G. Actual Errors 

Finally, for those particle orbits included in the 
highly accurate frequency-domain self-fo rce c alculations 
of Detweiler, Messaritaki, and Whiting jl69l | and Diaz- 
Rivera et at. [l82l |. I can compute the actual self- force 
errors as 
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(actual) 
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self - 
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(38) 



Here I'm implicitly taking the published results as "ex- 
act", i.e., I'm assuming that they're computed much 
more accurately than my computations. This is true for 
most, though not all, of the numerical results presented 
in this paper. In particular, for the test case considered in 
scction|Vl Detweiler, Messaritaki, and Whiting [l69| give 
the self force as F^df = 1.378 448 28 x IQ-^q^/M, with 
an estimated uncertainty of Admw = 2 x \Q~'^^q^ /M 
(0.015 ppm). I consider any "actual errors" defined 
by ([55]) which are less than 3Admw to be unreliable. 



H. Summary 

In summary, the numerical computation of the self 
force involves the following steps: 

1. Numerically solve the wave equation ^ for a suit- 
able set of m), using either double or long-double 
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floating-point arithmetic. (All subsequent steps use 
double floating-point arithmetic.) 



2. Calculate the regularized self-force modes -Ff,reg 
and their internal error estimates S F^^^^^'^^^^'^ for the 
corresponding set of £, using (fTS]) , (fTBI) , and ([32|) . 



ternal error estimate (^i^golf "JJj^^' 
one of the definitions (l33l). 



using (j22bp and 



4. Perform the tail fit to determine the coeffi- 
cients {cp}, then calculate the tail force Fscif,taii and 
its internal error estimate <5-Fsc[f tar^'' using 
(PU]). and one of the definitions ([M)) . 



5. Compute the self force -Fsoif and its internal error 



estimate 6F^ 



(internal) 



initions (|35p . 



self 



using p2ap and one of the def- 



V. NUMERICAL RESULTS 
A. Numerical Parameters 

As a test case I take Vp — lOM; the particle's orbital 
period is 2tt/u!p sa 199 M. In most cases I compute the 
regularized self-force modes ^f,rcg and their internal er- 
ror estimates SF^^^^^'^'^^^^ (step [2] in the summary of sec- 
tion llVH)) for < 30, ^ = 35, and £ = 40. Taking into 
account that aim — (and thus the wave equation ^ is 
trivial) if £ — m is odd, this set of £ gives a total of 295 
distinct (£, m) for which the wave equation ^ must be 
solved numerically (step[T|). 

To test the numerical computation over a wide range 
of cost /accuracy tradeoffs, I perform steps [T] and [2] for 
each of the combinations of the numerical-accuracy pa- 
rameters (floating-point precision, AMR error tolerance, 
and FMR refinement factor) shown with a "-y/" symbol 
in table El 

For a given AMR error tolerance euo, the regular- 
ized self-force modes ^i,rog have internal error esti- 
mates (5F^''^"*g which vary over almost 4 orders of mag- 
nitude over the range of £ I compute (this can be seen in 
figure [S]). This suggests that a better ratio of accuracy 
to computational cost might be obtained by applying a 
larger amount of FMR to those modes with the largest in- 
ternal error estimates, and a smaller amount of FMR (or 
none at all) for those modes with relatively small internal 
error estimates. For the long-double eite = 10~^^ calcu- 
lation, I thus also perform a further "playbaek23" calcu- 
lation which uses FMR by a factor of 3 for 12 < £ < 15, 
22 < £ <30, £ ^ 35, and ^ = 40 (these are shown with a 
"(V)" symbol in table El, and FMR by a factor of 2 for 
the other £. 

The required problem domain size for the numerical 
solution of the wave equation (jH]) is set by how long 



it takes the incorrect-initial-data perturbation to decay 
below the numerical error level. This size needs to 
be larger for smaller £ (where the perturbation decays 
more slowly) and for greater accuracy (smaller AMR 
error tolerance euo and/or larger FMR refinement fac- 
tor). For example, figure [3] shows the time dependence of 

-^0 reg (where the perturbation decays very slowly) and 

of -fi^rog (where the perturbation decays fairly rapidly); 
this latter case is is qualitatively similar to those of the 
other -P/^gg with £ > 0. Based on trial experiments 
with different problem-domain sizes, I have adopted the 



problem-domain sizes given in table ED 

Because of the very slow decay of £ — perturbations 
in Schwarzschild spacetime (this is visible in figure [3]) , I 
use very large problem-domain sizes for £ ~ to ensure 
that fI^^.I can be computed very accurately. Since I use 
an AMR numerical scheme jl78l | where the numerical 
evolution's computational cost is strongly concentrated 
near the particle worldline, even the very large £ ~ 
problem-domain sizes still only contribute a small frac- 



Numerical- Accuracy Parameters 



Floating-Point 



FMR Refinement Factor 



Precision 


£lte 


1 


2 


3 4 6 8 


double 


10-" 


V 


V 




double 


10-^5 


V 


V 




double 


10-16 


V 


V 


■ Lj_h 


long-double 


10-14 


V 


V 




long-double 


10-^5 


V 


V 




long-double 


10-16 


V 


V 




long-double 


10-1^ 


V 


V 




long-double 


10-18 


V 


V 




long-double 


10-19 


V 


V 


(%/) 



^ = < £ < 30, ^ = 35, ^ = 40 

V = < ^ < 30, ^ = 35, ^ = 40; 

record-playback error estimate can be computed 

(^) = 12 < ^ < 15, 22 < ^ < 30, £ = 35, ^ = 40 

^1 = solution of the wave equation ((9]) is seriously 
affected by floating-point roundoff errors 

TABLE V. This table lists the numerical-accuracy parame- 
ters (the floating-point precision, AMR error tolerance eito, 
and FMR refinement factor), and £ for which I have nu- 
merically solved the wave equation ^ (step [1] in the sum- 
mary of section [IV H|) and computed the regularized self- force 
modes -Ff,rcg (step [21). (This latter computation always uses 
double floating-point precision.) The shaded cells mark pa- 
rameters where the solution of the wave equation Q is seri- 
ously affected by floating-point roundoff errors. (Results from 
these parameters are plotted as the "double (bad)" points in 
figure [HI ) 
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tion (typically < 15%) of the total computational cost of 
the self-force calculatioiiPlcJ 



Vertical Scale 
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FIG. 3. [Color online] This figure shows the decay of the 
regularized self-force modes ^o.^ig (top) and ^"itJ^reg (bottom) 
towards their late-time values, for the £ito = 10~^^ "record" 
evolution. The vertical scale changes at each vertical dashed 
line, zooming in from left to right by the factors shown above 
each plot. (The changing vertical scale also accounts for the 



apparent increase in the difference Ff^J^ — F, ^' at later times; 
this difference is actually almost time-independent.) In the 



?(-) 



0,reg 



plot the horizontal scale changes from logarithmic for 
ischw < 30 OOOA/ (shown as the upper of the two rows of time 
labels below the plot) to linear for tschw > 30 DOOM (shown 
as the lower of the two rows of time labels below the plot). 



The problem-domain sizes shown in table I VII suffice 
to ensure that each Fi^^eg is time-independent to well 
within the numerical errors by the end of its numerical 
evolution. All the results reported here use Pj:"^^^ values 
sampled lOM before the end of the evolution. 

To further explore cost / accuracy tradeoffs in self- force 
calculations, for each combination of numerical-accuracy 
parameters given in tablc|V]l have computed the numeri- 
cal force i^scif,num (step[3]in the summary of section lTV Hp . 
performed the tail fit and computed the tail force i^scif.taii 
(step [4]), and computed the self- force (step [5]) for each for 
the numerical-force and tail-fit parameters shown in ta- 
bic IVlTl 



B. Overview of the Numerical Results 

Figure 3] shows the regularized self-force modes Ff^rog 
for the most accurate "record" evolution {e\tc = 10"-^^). 
Notice that the large-^ modes are very closely approxi- 
mated by the tail fit ([23|) : I discuss this further in sec- 
tion |VDl 

Figure [5] shows the regularized self-force modes' in- 
ternal error estimates iJi^/'^^'g"^"^'^ for the record, play- 
back2, and playback23 £itc = 10~^^ evolutions. As £ 
increases, the solutions of the wave equation ^ oscillate 
more rapidly in space and time, so the finite-difference 
truncation errors for any fixed numerical resolution in- 
crease rapidly. The AMR algorithm responds to this by 
decreasing Avu^nin in discrete factor-of-2 steps, each of 
which decreases the global finite-difference truncation er- 
ror of the solution by very close to a factor of 2* = 16. 
This accounts for the "stepped" appearance of the er- 
ror estimates visible in figure El (Notice also that as 
intended, the playback23 error estimates show much less 
dynamic range than the record or playback2 estimates.) 

Problem-Domain Sizes 
£ AMR Error Tolerance £ito Problem-Domain Size (Af) 





1 
2 

3-4 
5-30 
35 
40 



10 
10 



-14 
17 



10" 
10" 



10" 
10" 



30 000 
100 000 
5 000 
1000 
500 
400 
400 
400 



With an AMR scheme of this type, the total cost of a numer- 
ical evolution at a given accuracy grows only linearly with the 
problem-domain size, rather than quadratically as would be the 
case for a unigrid scheme. 

For an AMR scheme such as mine there would be little benefit in 
using a timelike inn er bo undary (and boundary condition) of the 
type used by Haas [l76j |: while this could remove almost half of 



TABLE VI. This table shows the problem-domain size used 
for each numerical evolution. 



the problem domain, the region removed would be distant from 
the particle worldlinc, so its removal would only save a small part 
of the total computational cost of the numerical evolution. 
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Depending on £, the AMR algorithm uses between 5 
and 8 levels of 2:1 mesh refinement for these evolutions. 
For the record evolution, the grid resolution of the finest 
refinement level varies from A^/128 to M/2048 depending 
on £] the resolutions for the playback2 and playback23 
evolutions are correspondingly finer. 

The speedup factor of the AMR algorithm over 
an equivalent-resolution (and thus roughly equivalent- 
accuracy) FMR evolution is typically 30 to 40 for evo- 
lutions using an £ = problem-domain size of 30 OOOAf , 
and 200 to 400 for evolutions (such as the eite = 10~^^ 
ones) using an £ = problem-domain size of 100 OOOA/. 



C. Validation of the Error Estimates 



sion are generally generally consistent with the record- 
playback error estimates Sf'^'^^ , and are somewhat con- 
servative (the internal error estimates tend to overesti- 
mate the record-playback error estimates). 

For evolutions done in double floating-point precision 
the internal error estimates are similarly consistent and 
conservative for the "good" parameters which are not 
shown as shaded in table |V] However, for the "bad" pa- 
rameters which are shown as shaded in table |V] the in- 
ternal error estimates scatter widely about the record- 
playback error estimates, often deviating by up to two or- 
ders of magnitude. This is due to fioating-point rounding 
errors contaminating the numerical solution of the wave 
equation ((9)) (step [T] in the summary of section IIV H|) . 



In this section I present numerical tests to validate 
the internal error estimates described in section IIV El 
against the record-playback error estimates described in 
section HV Fl and (in those cases where the actual errors 
are known) to validate the record-playback error esti- 
mates against the actual errors. For the comparisons of 
internal with record-playback error estimates I use the 
results from all of the numerical- accuracy parameters for 
which a record-playback error estimate can be computed 
(these parameters are shown in table |V|. To prevent in- 
accuracies in the tail fit from contaminating the error 
estimates, in this section I consider only the highest- 
accuracy set of numerical-force and tail-fit parameters 
shown in table IVIII (i.e., those in in the last row of the 
table): the numerical force sums modes up to X = 30, 
the tail fit includes the modes 20 < £ < 30, £ = 35, and 
£ = 40, and the tail fit fits either {c4, cg} (with C2 given 
analytically by (HU) or {c2,C4,C6}. 



Individual Regularized Self-Force Modes Fi^. 



Figure [6] shows a scatterplot of the the record-playback 



error estimates 5F, 



(r-p) 



versus the internal error esti- 



mates 



The internal error estimates 5F, 



(internal) 



for evolutions done in long-double floating-point preci- 



Numerical-Force and Tail-Fit Parameters 



Analytical C2 



Tail-Fitted C2 



2. Self-Force Numerical Sum Fseif,, 



Figure [7] shows a scatterplot of the numerical-force 



record-playback error estimate SF^ 



versus the in- 



ternal error estimate SF^ 



(internal) 



SGlf,num 

soif.num ' ^^'^ latter computed 
using each of the definitions ([55]) . The arithmetic- 
sum internal error estimate p3bp systematically overes- 
timates the record-playback error estimate by a factor 
of '--^ 3.5. In contrast, the quadrature-sum internal er- 
ror estimate (|33ap is quite accurate^ Based on this, I 
adopt the quadrature-sum internal error estimate p3ap 
hereinafter. 
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K {£} in tail fit {4} {4, 6} {2, 4} {2, 4, 6} FIG. 4. [Color online] This figure shows the regularized self- 
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10-15 
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20-30, 35, 40 



V 
V 
V 
V 



V 
V 
V 
V 



r ^ j-ifintcrnal) 

force modes -f/rcg 
then logarithmic from £ 



The £ scale is linear from £ = to 3, 
3 to 40. 



TABLE VII. This table shows the numerical-force and tail-fit 
parameters. K is the maximum £ included in the numerical 
force. 



^® Notice that the set of modes considered here includes the "double 
(bad)" regularized self-force modes plotted in figure [6] and dis- 
cussed in section FV C 1! Evidently the averaging inherent in sum- 
ming 33 of these modes greatly reduces the effects of the floating- 
point roundoff error contamination of the individual modes. 
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3. Self -Force Tail Sum Fseif.taU 

Figure [8] shows a scatterplot of the tail-force record- 
playback error estimate ^F^l\f\^{\ versus the internal er- 
ror estimate {the latter computed using each 
combination of the definitions p4|) and fitting either 
{c4,C6} or {c2,C4,C6}). When fitting {c4,C6} (with C2 
given analytically by (US])), for both double and long- 
double floating-point precision the worst-case-of-3*^-trials 
internal error estimate (|34bp tends to systematically over- 
estimate the record-playback error estimate, while the 
statistical internal error estimate (|34ap is much more ac- 
curate. Based on this, I adopt the statistical internal 
error estimate (|34ap hereinafter. For both precisions, the 
internal error estimates change from being overestimates 
to underestimates of the record-playback error estimates 
for the two smallest-error points. 

When fitting {02,0^,00}, the long-double internal er- 
ror estimates are still consistent and somewhat conser- 
vative, but the three smallest-error internal error esti- 
mates scatter widely about the corresponding record- 
playback error estimates. This appears to be due to 
the ill-conditioning of the {c2, C4, cq} tail fit (cf. discus- 
sion in section lTV Dl particularly table lIV[) amplifying the 
floating-point rounding errors in the individual regular- 
ized self-force modes F(,^cg- 

I discuss further "quality control" checks based on the 
tail fits' a-nd residuals in section IVDl 



terplot of the actual self-force errors SF^^^f versus the 
self-force internal error estimates (^^j^p latter 

computed using each combination of the definitions psp ) 
and for tail fits to {04, ce} or {c2, C4, cg}). All of the error 
estimates are fairly accurate. 

5. Record- Playback Error Estimates 

The self- force actual errors can also be used to validate 
the record-playback error estimates. Figure [10] shows a 
scatterplot of the actual self-force errors versus 

the record-playback error estimates SF^^^^\ The actual 
errors are very similar to the record-playback error esti- 
mates. 



D. The Tail Fits 

The quality of a tail fit can be assessed via the fit's x^- 
if lies outside the "plausible" range [X2 5%' ^97 5%]' 
where is the /3% percentile of the x^ distribution for 
the appropriate number of degrees of freedom (here 11 for 
fitting {c4, cg}, or 10 for fitting {c2, C4, ce}), we can reject 
the null hypothesis that the the nonzero fit residuals are 
solely due to (independent Gaussian-distributed) random 
errors of magnitudes given by the individual modes' in- 
ternal error estimates ""^'^ A corollary of this is 



4- Overall Self-Force 



As discussed in section IIV G[ I can compute the ac- 
tual error of the overall self-force by comparing my cal- 
culations against previously published highly-accurate 
frequency-domain calculations. Figure [H] shows a scat- 
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FIG. 5. [Color online] This figure shows the regularized self- 



force modes' internal error estimates SF1\ 



{internal) 



£,rcg 



for the eito 



10 evolutions. The £ scale is linear from = to 3, then 
logarithmic from £ = 3 to 40. 
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FIG. 6. [Color online] This figure shows a scatterplot of the 
record-playback error estimate SF^'^'^^ versus the internal er- 
ror estimate S F^^^^^'^'^'^^K The solid and dashed lines show the 
cases where the record-playback error estimate is identical 
to or an order of magnitude larger/smaller than the inter- 
nal error estimate, respectively. The relative-error scales are 
relative to the overall self force -Fseif- 
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that the statistical tail-fit error estimate (|34aP becomes 
umeliable. 

Figure fTD shows the tail-fit for fitting both {04, cg} 
and {c2, C4, cg} for each set of numerical-accuracy param- 
eters listed in table |Vl For effective error tolerances 
gite.e ff ^. 10^^° the fits are very well-behaved: is 

Aft rpsr and 

decrease oc eft^'^^ff 



smallFI and both the RMS residuals 



the self-force actual errors Jft 



(actual) 



self 

Eite.eff decreases, as expected for a characteristic evolution 
scheme with 4th/6th order global/local finite differencing 
accuracy. 

However, for £itc,cff ^ 10^^° the fits show several un- 
desirable characteristics: increases, the RMS residu- 
als either begin to increase with decreasing eito,cff (dou- 
ble floating-point precision) or decrease at a slower rate 

than cx eft{,"^(,ff (long-double floating-point precision), and 
the actual errors cither increase for the very small- 
est eito,Gff (double floating-point precision), level off as 
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FIG. 7. [Color online] This figure shows a scatterplot of the 
numerical-force record-playback estimate ^^'a'ei/'nmn versus the 

internal error estimate '5J"gciTn™''' (^^^ latter computed using 
each of the definitions (|33|) ). The solid and dashed lines show 
the cases where the record-playback error estimate is identical 
to or an order of magnitude larger/smaller than the internal 
error estimate, respectively. The relative-error scales are rel- 
ative to the overall self force Fsdf. 



In fact, X2 5%- This is due to the individual regularized 

self-force modes' internal error estimates ^pi'^*-'^'^'^^^) systemati- 
cally overrestimating the actual numerical errors (cf. figure[6]and 
discussion in section [V C 11 1. 
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FIG. 8. [Color online] This figure shows a scatterplot of 
the tail-force record-playback estimate 5^jei/taii versus the 
internal error estimate ^-Fjelt'tan"'' i^^^ latter computed using 
each of the definitions (|34p 'l. for tail fits to {c4, ce} (top) and 
{c2,C4,C6} (bottom). The solid and dashed lines show the 
cases where the record-playback error estimate is identical 
to or an order of magnitude larger/smaller than the inter- 
nal error estimate, respectively. The relative-error scales are 
relative to the overall self force -FLeif. 
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FIG. 9. [Color online] This figure shows a scatterplot of the 
self-force actual errors Sp'^^'j}^^^^ versus the internal error esti- 

sclr 

mates j^J^'J^*'"'"^') (•^jjg latter computed using each of the defi- 
nitions (|35[l ). for tail fits to {c4, cg} (top) and {c2, C4, cg} (bot- 
tom). The solid and dashed lines show the cases where 
the actual errors are identical to or an order of magnitude 
larger/smaller than the internal error estimates, respectively. 
The actual-error values are unreliable in the shaded region of 
each plot. The relative-error scales are relative to the overall 
self force -Fscif. 



^itc.cff decreases (long-double floating-point precision, tail 
fit to {ci,Ce} only), or decrease at a slower rate than 
cx eft{^eff (long-double floating-point precision, tail fit to 
{c2, C4, cg}. These effects are due to floating-point round- 
off errors contaminating the various steps of the calcula- 
tion. 



E. Cost/ Accuracy Tradeoffs 

There are a number of cost /accuracy tradeoffs inherent 
in the choice of the various numerical parameters in the 
self-force computation. 

The computational cost is overwhelmingly dominated 
by the numerical solution of the wave equation (jH]) (stcp[T] 
in the summary of section |IVH[) . and is determined by 
the combination of the set of £ for which the wave equa- 
tion ([9]) is solved, and the problem-domain sizes, floating- 
point precision, AMR error tolerances eue, and FMR re- 
finement factors (if any) used in that solution. 

In general, the problem-domain size should be chosen 
just large enough to render the errors induced by the re- 
maining time dependence of ^£,reg small in comparison 
to other numerical errors. The problem-domain size re- 
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FIG. 10. [Color online] This figure shows a scatterplot of the 
self-force actual errors gpj-^'^'-^^^^ 
error estimates 5F_ 



self versus the record-playback 
The solid and dashed lines show the 



cases where the actual errors are identical to or an order of 
magnitude larger/smaller than the record-playback error esti- 
mates, respectively. The actual-error values are unreliable in 
the shaded region of each plot. The relative-error scales are 
relative to the overall self force Fseif- 
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quired to ensure this varies with the magnitude of the 
other numerical errors, being larger for higher accura- 
cies (smaller errors). The required problem-domain size 
also varies strongly with being much larger for small (. 
(cf . discussion in section IV Ap . For simplicity, in this 
work I have only adjusted the problem-domain sizes at 
the very coarse level shown in table IVIl A more careful 
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adjustment would substantially improve the efficiency of 
the computation. 

For my computational scheme, AMR and FMR are 
of almost equal efficiencies. That is, the grid structure, 
computational cost, and accuracy attained from an evo- 
lution using AMR with error tolerance eitc and FMR by 
a refinement factor of iVtmr arc all very similar to those 
obtained from a purely-AMR evolution using the corre- 
sponding effective error tolerance eito,off = £ito/A^fmr- 
suming that the problem-domain sizes are large enough 
so that the remaining time dependence of ^£,reg isn't a 
significant contributor to the overall error budget, the pa- 
rameter space for cost-performance tradeoffs in numeri- 
cally solving the wave equation ^ can thus be simplified 
to just the effective error tolerance eito,off- 

The accuracy of a self-force computation is then de- 
termined by the combination of the set of £ for which 
the wave equation ^ is solved, the effective error toler- 
ance eito,off of this solution, K (the maximum £ included 
in the numerical force), the set of £ used in fitting the 
tail series (|23p , and the set of orders p and corresponding 
coefficients {cp} included in this series. This is a large 
parameter space; for present purposes I restrict consider- 
ation to those parameter combinations listed in tables IVl 
andEIIl 

For present purposes, it's useful to quantify the com- 
putational cost of an evolution by the total number of 
diamond cells integrated by the AMR algorithm. This is 
closely proportional to the overall CPU time used, with 
the constant of proportionality (the CPU time per di- 
amond cell) being about 1.5 (3.0) microseconds per di- 
amond cell for double (long-double) fioating-point pre- 
cision on the processors used here. Figure [T^] gives an 
overview of the cost-accuracy tradeoffs for the highest- 
accuracy set of numerical-force and tail-fit parameters 
shown in table IVIII (i.e., those in in the last row of the 
table): the numerical force sums modes up to ii' = 30, 
the tail fit includes the modes 20 < £ < 30, £ = 35, and 
£ = 40, and the tail fit fits either {c4, cg} (with C2 given 
analytically by p5|) ') or {02,04, cq}. It should be noted 
that the costs shown in this figure are for computations 
with very conservative (large) problem-domain sizes, and 
the wave equation © solve for a large set of £. The costs 
could be greatly reduced with only a minor impact on 
the accuracy by using smaller problem-domain sizes and 
a smaller set of £. 



FIG. 11. [Color online] This figure shows the (upper plot), 
RMS residuals || AF£,regj|j.j^j, (middle plot), and self-force ac- 
tual error gp^^^^^^^'i (lower plot) for each set of numerical- 
accuracy parameters listed in table [V] In the plot, the two 
pairs of dashed lines show the 2.5% and 97.5% percentiles of 
the distribution for 11 degrees of freedom (appropriate for 
fitting {c4,C6}) and for 10 degrees of freedom (appropriate 
for fitting {c2, C4, cg}). In the RMS-residual plot the diagonal 
line shows the s^^ ^^ff scalings expected for my characteristic 
AMR algorithm [17g|. In the actual-error plot, the actual- 
error values are unreliable in the shaded region. 



Results for the Self-Force 



Tables IVIIII and IIXI give the main results of my com- 
putations for the self-force and its error estimates. 

These results are fully consistent at the 0.1 ppm level 
with the highly accurate frequency-domain results of De- 
tweiler, Messaritaki, and Whiting |l69l |. Below O.I ppm 
my results' error estimates become increasingly unre- 
liable due to fioating-point roundoff errors, and below 
0.045 ppm the finite accuracy of the Detweiler, Messar- 
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itaki, and Whiting jl69| results (their quoted error esti- 
mate is 0.015 ppm) begins to affect comparisons with my 
results. 



VI. CONCLUSIONS 

This work demonstrates that the use of character- 
istic AMR can dramatically improve the efficiency of 
time-domain self-force calculations using the Barack-Ori 
mode-sum regularization formalism!^ 

I find that the tail-fit basis {fp} is very ill-conditioned 
if many terms in the tail series (|23|) arc fit simultaneously. 
Fortunately, normalizing the basis functions to have sim- 
ilar magnitudes mostly alleviates this ill-conditioning. 

Past self-force calculations have often used "record- 
playback" error estimates derived from comparing a pair 
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FIG. 12. [Color online] This figure shows an overview of 
the cost-accuracy tradeoffs for the highest-accuracy set of 
numerical-force and tail-fit parameters shown in table IVIII 
Notice that there are two different CPU-time scales: the up- 
per (outer) scale is for evolutions in long-double floating-point 
precision, while the lower (inner) scale is for evolutions in dou- 
ble floating-point precision. The solid line shows the expected 
scaling ^Fjj'if""') ociV-2. 



of different-resolution calculations. Here I present, and 
validate as quite reliable, a set of internal error estima- 
tors which can be used within a single self-force calcula- 
tion (whether AMR or unigrid) to estimate the accuracy 
of individual regularized self- force modes -Fi.rog, and the 
numerical force, tail force, and overall self-force derived 
from them. 

In their pioneering calculation of the gravitational self- 
force acting on a mass particle on a circular geodesi c or- 
bit in Schwarzschild spacetimc, Barack and Sago |l77| 
use the arithmetic-sum formula (j33b| to combine the 
numerical-force error given the (record-playback) error 
estimates of the individual modes. Here I show that (at 
least for my results) this formula is unnecessarily conser- 
vative, and that the quadrature-sum formula pSap pro- 
vides a better approximation to the numerical-force error 
over a wide range of overall computational acc urac i es. 

Like other researchers (sec, for example, |l62l . Il77l | 
and references therein, and the references cited in foot- 
note I13p , I find excellent agreement between time- and 
frequency-domain calculations of the self-force. Here I 
demonstrate this agreement down to the 0.1 ppm accu- 
racy level. Because the time- and frequency-domain cal- 
culations are structured so differently, this high-precision 
verification of their agreement provides a strong confir- 
mation of the correctness of both calculations, and im- 
plicitly of their respective theoretical formalisms as well. 

The present work could (should) be extended in sev- 
eral directions. Apart from the technical limits of the 
relatively coarse adjustment of the problem-domain size, 
one obvious extension would be to consider the electro- 
magnetic and/or gravitational self-force. Another possi- 
bility would be to generalize the current finite differenc- 
ing scheme to handle non-circular particle orbits. This 
would be straightforward, albeit somewhat te dious , us- 
ing techniques suc h as those described by Haas |l76t and 
Barack and Sago |l77l |. Once non-circular orbits can be 
handled, it should then be possible to move to full orbit- 
correction calculations of the type suggested by Gralla 
and Wald |l57l section 7] . This will be very computation- 
ally demanding (it might benefit from further increasing 
the order of the finite differencing), but should provide 
valuable information about 0(^^) radiation-reaction ef- 
fects. 

The generalization of this work to particle orbits in 
Kerr spacetime would also be very valuable, but would 
demand a major reorganization of the mathematical and 
computational structure. 
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Tail fit fits coefficients {c^, cq} (c2 is given analytically by the circular-orbit formula ( |25l l): 
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TABLE VIII. This table shows the main results of the self-force calculations for the case where the tail fit fits only {c4, ce}. For 
each calculation, the table shows the AMR error tolerance eitc used in numerically solving the wave equation (j9}, whether this 
numerical solution is purely AMR ("record") or also uses FMR ( "playbackA''" for some TV), the computed self- force -Fscif, for 
the tail fit, the internal error estimate 5F^^^^'"^'^'^^\ the record- playback error estimate 5F^^^^\ and the actual error ^^J^^"^^*"^') 
Each error estimate or error is shown both as an absolute value, and as a relative value in parts per million (ppm) relative to the 
overall self- force. The shaded rows have very large tail-fit x^, so their internal estimates may be unreliable. For compariso n, th e 
last row of this table shows the highly accurate frequency-domain value calculated by Detweiler, Messaritaki, and Whiting [l69l |. 
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TABLE IX. This table shows the main results of the self- force calculations for the case where the tail fit fits {c2,C4,C6}, as 
might be the case for a non-circular particle orbit. For each calculation, the table shows the AMR error tolerance eue used 
in numerically solving the wave equation @, whether this numerical solution is purely AMR ("record") or also uses FMR 
( "playbackA'^" for some A''), the computed self-force Fsdf, for the tail fit, the internal error estimate (J-Fj^JJf , the record- 
playback error estimate SF^I^^^\ and the actual error 

^^Jactual)^ Each error estimate or error is shown both as an absolute 
value, and as a relative value in parts per million (ppm) relative to the overall self-force. The shaded rows have very large 
tail-fit so their internal estimates may be unreliable. For comparison, t he la st row of this table shows the highly accurate 
frequency-domain value calculated by Detweiler, Messaritaki, and Whiting [169| |. 
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